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GAUSSIAN  RAY  BUNDLES  FOR  MODELING 
HIGH-FREQUENCY  PROPAGATION  LOSS 
UNDER  SHALLOW-WATER  CONDITIONS 


1.  INTRODUCTION 


This  report  introduces  a  range-dependent  propagation  loss  model  designed  to  support  the 
simulation  of  high-frequency  reverberation  in  shallow  oceans.  For  simplicity,  this  discussion  is 
confined  ioNxlD  environments.  That  is,  environmental  quantities  such  as  the  sound  speed  and 
bottom  depth  vary  along  vertical  (r,  z)  planes  through  the  source  and  field  points,  and  out-of¬ 
plane  scattering  is  neglected.  The  model  is  based  on  Gaussian  ray  bundles  of  the  form 


0) 


where  Po  depends  only  on  the  source,  /’includes  losses  due  to  volume  attenuation  and  boundary 
reflections,  Zy  is  the  depth  along  a  central  ray,  pr  is  the  horizontal  slowness,  and  cr,  an  effective 
standard  deviation  or  half-beamwidth,  is  given  by 


cr  =  ^max(Az,8/l) . 


(2) 


Here,  Az  is  the  change  in  ray  depth  at  constant  range  due  to  a  change  in  source  angle  A  61),  and  X 
is  the  wavelength  at  the  field  point. 


If  Az  >  SX,  (T  confines  most  of  the  energy  to  the  region  between  adjacent  rays,  while  the  SX 
minimum  aperture  prevents  the  acoustic  pressure  from  growing  too  large  near  caustics,  where 


lim  =  0. 

A0o->O  A^q 


(3) 


More  sophisticated  caustic  corrections  are  often  difficult  to  apply  to  realistic  ocean 
conditions  due  to  the  numerical  instability  of  various  parameters  or  the  violation  of  basic 
assumptions.  For  example,  Ludwig's  (reference  1)  uniform  asymptotic  expansion  involving  the 
Airy  function  Ai  and  its  derivative  Ai'  applies  to  smooth  caustics.  There  is  a  smooth  transition 
between  the  oscillatory  and  damped  behavior  with  large  but  finite  values  on  the  caustic  itself  The 
expansion  requires  equation  (3)  to  have  a  well-behaved  first-order  zero.  This  assumption  is 
violated  along  rays  that  graze  an  ocean  boundary  because  the  limit  in  equation  (3)  fails  to  exist. 

Equation  (1)  is  based  on  conservation  laws  and  appears  to  be  consistent  with  Fresnel 
volume  ray  tracing  (reference  2).  Equation  (2),  on  the  other  hand,  was  discovered  by  accident 
during  the  development  of  a  Gaussian  beam  model. 
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A  Gaussian  beam  (reference  3)  is  defined  along  a  central  ray  by 


P(s,  rj)  =  Ayjc(s)/[rq(s)]  x  exp(-/®  {7(5)  +  0.5[/?(5)  /  q{s)\Tf  })  .  (4) 

Again,  r  is  the  horizontal  range,  s  is  the  arc  length  along  the  ray,  77  is  the  normal  distance  from  the 
central  ray,  ^  is  an  arbitrary  constant,  o  is  the  radian  frequency,  c(5)  is  the  sound  speed,  and  p{s) 
and  q{s)  are  the  beam  curvature  and  beamwidth,  respectively.  The  curvature  and  beamwidth 
satisfy  a  pair  of  ordinary  differential  equations  that  are  integrated  with  the  standard  ray-tracing 
equations. 

Gaussian  ray  bundles  are  somewhat  easier  to  determine  than  Gaussian  beams  since  <j  is 
determined  from  adjacent  rays  at  the  field  point,  not  by  integrating  differential  equations.  It  is  not 
known  why  the  8A  in  equation  (2)  performs  as  well  as  it  does.  Hence,  the  bundles  should  be  used 
with  caution  until  cr  can  be  placed  on  a  firmer  foundation,  or  until  another  model  can  be  shown  to 
provide  a  practical  solution  at  the  high  frequencies  of  interest. 

Several  aspects  of  Gaussian  ray  bundles  resemble  the  numerical  integration  of  the  one-way 
wide-angle  Maslov-Chapman  wavefield  representation  (reference  4).  For  example,  transforming 
the  Helmholtz  equation  from  depth  dependence  to  vertical  slowness  leads  to  a  modified  eikonal 
equation.  The  solution,  a  Legendre  transform  from  depth  to  vertical  slowness,  serves  as  an 
algorithm,  equation  (3 1),  to  extrapolate  ray  bundle  travel  time  to  a  field  point.  An  important 
difference  is  that  the  present  numerical  integration  of  the  pressure  amplitude  is  based  on  power 
addition,  not  on  coherent  transforms. 

Hardy  (reference  5)  suggested  that  a  Gaussian  form  involving  an  experimentally  determined 
standard  deviation  could  be  used  to  treat  caustics.  Although  critics  argued  that  one  would  have 
to  conduct  a  separate  experiment  for  every  case  of  interest,  it  now  appears  that  the  single 
expression,  equation  (2),  is  sufficient.  This  will  be  demonstrated  by  several  heuristic  arguments  in 
section  2. 

Section  3  contains  two  propagation  loss  test  cases.  The  first  compares  results  with  those  of 
the  Navy  standards  PE  version  3.4  (reference  6)  and  ASTRAL  version  4.2  (reference  7)  at  25  and 
10  kHz  for  a  classic  convergence  zone.  Then,  the  more  academic  models  EFEPE  (reference  8) 
and  OASES  (reference  9)  are  used  to  provide  propagation  loss  at  1000  Hz  for  a  shallow- water 
environment.  If  Gaussian  ray  bundles  compare  well  with  these  “standards”  at  the  lower 
frequencies,  they  should  perform  well  at  the  higher  frequencies,  since  the  increase  in  frequency 
only  improves  the  high-frequency  approximation  invoked  by  ray  theory. 

The  ray  segments  used  by  the  classical  ray-tracing  section  are  based  on  triangular  sectors  in 
which  the  inverse  sound  speed  squared  is  a  linear  function  of  range  and  depth.  The  specific 
implementation,  being  somewhat  different  than  previous  implementations,  is  described  in  the 
appendix. 
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2.  GAUSSIAN  RAY  BUNDLES 


The  original  intent  of  this  work  was  to  evolve  an  iV  x  2Z)  ray-tracing  program  into  a 
Gaussian  beam  model.  Development  by  evolution  often  gives  additional  insight.  In  this  instance, 
it  led  to  Gaussian  ray  bundles.  Since  rigorous  proofs  are  not  currently  available,  the  justification 
relies  on  heuristic  arguments.  This  section  begins  with  a  brief  discussion  of  classical  ray  theory  as 
a  vehicle  to  introduce  the  notation.  Gaussian  ray  bundles  are  constructed  from  classical  rays  in  a 
number  of  ways,  four  of  which  are  described  below.  The  basic  technique  applies  to  a  benign 
environment  in  which  there  are  sufficient  classical  rays  of  a  desired  ray  type  to  form  well-behaved 
differences.  The  second  technique  is  di  fallback  procedure  in  case  there  is  only  one  classical  ray  of 
a  particular  type  in  the  region  of  interest.  Virtual  rays  are  defined  by  unfolding  classical  rays  to 
model  propagation  in  the  vicinity  of  the  ocean  boundaries.  The  last  technique  modifies  the  basic 
formula  to  model  propagation  at  short  ranges. 

First  consider  the  fan  of  n  rays  in  figure  1 .  Parameters  that  describe  the  v-th  ray  are 
subscripted  by  ( ),^  while  ( )o  refers  to  a  point  source  with  radian  frequency  (o  at  (0,  zo). 
Accordingly,  Po  denotes  the  acoustic  pressure  at  a  reference  distance  ro  from  the  source.  The  v-th 
ray  is  launched  with  source  angle  Ovst  and  crosses  the  field-point  range  r  at  depth  Zv  and  travel 
time  Tv. 

The  horizontal  slowness  and  vertical  slowness, 


P..v  =  cos  ejc. 


Figure  1.  Fan  of  Acoustic  Rays  Launched  from  a  Point  Source 
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are  defined  in  terms  of  the  horizontal  inclination  angle  6v  and  the  sound  speed  Cy.  Losses  due  to 
volume  attenuation  and  boundary  reflections  are  simulated  by  the  pressure  ratio  Fy  and  phase  shift 
0v.  Recall  that  the  vertex  velocity  is  a  ray  invariant  in  range-independent  environments. 


Typical  Nx  2D  ray-tracing  programs  require  two  types  of  spatial  interpolations  to  model 
realistic  ocean  conditions.  The  first  provides  environmental  data  at  discrete  ranges  along  a 
vertical  plane.  The  second  interpolation  provides  environmental  data  at  arbitrary  points  in  the 
vertical  plane  so  that  ray  trajectories  can  be  evaluated  analytically.  In  particular,  if  the  inverse 
sound  speed  squared  is  found  by  linear  interpolation  within  triangular  sectors,  the  ray  trajectories 
become  parabolas.  The  implementation  used  here,  being  somewhat  different  than  others,  is 
described  in  the  appendix. 

The  numerical  implementation  of  geometrical  spreading  uses  finite  differences  to  estimate 
Azy,  the  increment  in  ray  depth  at  constant  range  due  to  an  increment  in  source  angle  AOy^.  For 
this  discussion,  it  is  sufficient  to  use  the  central  differences 

=(^v+i -■^v-i)/2,  (6) 

(7) 


for  v=  2,  3, ...,  « - 1,  and 

> 

A^,  0  =  02,0  ~  ^1,0  ’ 

A^n.O  ~  0n.O  ~  ^n-1,0 


(8) 

(9) 

(10) 

(11) 


at  the  endpoints. 

Later  examples  will  demonstrate  that  the  ocean  boundaries  will  make  the  above  definitions 
inappropriate  for  some  Azy.  In  anticipation  of  these  boundary-induced  complications,  two  depths 
and  z^^^  that  bound  Zy  are  now  introduced.  Also,  set 


Az  =  z*'^^  —  z 


rO) 


(12) 


Equation  (12)  will  be  equivalent  to  equations  (6),  (8),  and  (10)  if  the  z^J^  and  zf^  are  midpoints  of 
the  Zy.  That  is, 

=(Vi +^v)/2  (13) 
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for  v=  2,  3, n,  and 


(14) 

for  v=  1,  2, ...,«- 1.  The  remaining  depths, 

=(3^1 -•^2)/2  (15) 

and 

=(3^«  -  Vi)/2,  (16) 

are  obtained  by  reflecting  about  zi  and  z^’^  about  z„. 


Under  the  Nx2D  and  geometric-acoustic  approximations,  the  contribution  of  the  v-th  ray 
to  the  total  field  can  be  estimated  fi'om 


P  =r  PT 

^  V  O'*-  V 


p  rAz 

-4  r  V 


Pr,v,0^^  v,0 


■1/2 


exp(/<y7^  -t-zOv) , 


(17) 


if 


rO) 


<z  <z 


(2) 


or  if 


z 


(2) 


<z  <z 


(1) 


(18) 


(19) 


and  is  zero  otherwise.  Each  nontrivial  Py  corresponds  to  an  eigenray,  that  is,  a  ray  path  that 
connects  the  source  and  field  point.  The  underlying  concept  of  classical  ray  theory  is  that  the 
acoustic  pressure  P^  at  any  given  field  point  is  the  coherent  addition 

P.=ZP.  (20) 


of  eigenrays.  The  power  addition 

i^’/=zinr 

V 


(21) 


is  a  smoother  representation  and  better  suited  to  many  practical  applications.  Both  summations 
include  all  nonzero  values  of  Py 
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Figure  1  is  deceptively  simple  in  that  there  is  only  one  eigenray  per  field  point.  In  other 
words,  there  is  only  one  nonzero  Pv  In  general,  the  task  of  identifying  all  the  eigenrays  is 
difficult,  especially  in  a  range-dependent  environment.  To  demonstrate  this  fact,  it  is  convenient 
to  define  ray  type  as  the  number  of  surface  and  bottom  reflections  and  upper  and  lower  vertices. 
For  example,  the  primary  eigenrays  for  an  isovelocity,  flat-bottom  environment  consist  of  the 
direct  path,  the  one  surface  bounce  path,  the  one  bottom  bounce  path,  the  one  bottom  then 
surface  bounce  path,  and  the  one  surface  then  bottom  bounce  path,  illustrated  in  the  upper  portion 
of  figure  2.  The  upper  portion  of  figure  3  illustrates  the  continuous,  predictable  manner  in  which 
the  angles  of  these  ray  paths  vary  with  field-point  depth  at  a  0.463-km  (0.25-nmi)  range. 

The  possible  number  of  primary  eigenrays  in  a  range-dependent  environment  is  unlimited, 
and  the  task  of  identifying  all  the  eigenrays  can  be  formidable.  The  lower  portions  of  figures  2 
and  3  show  how  quickly  the  number  of  eigenrays  can  multiply  just  from  varying  the  bathymetry 
for  this  isovelocity  environment.  In  this  example,  the  direct  and  surface  bounce  paths  are 
unaffected,  but  the  number  of  bottom  bounce  paths  has  increased  fi'om  three  paths  to  nine  paths 
for  the  60.96-m  (200-ft)  field-point  depth  (figure  2).  The  arrival  angle  structure  (figure  3)  shows 
how  the  number  of  paths  changes  with  field-point  depth. 
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Figure  2.  Primary  Eigenrays  to  a  0.25-nmi,  100-ft  Field  Point  in  an  Isovelocity 
Environment:  (a)  Flat  and  (b)  Range-Dependent  Bathymetry 


ANGLE  (deg) 

-60  -40  -20  0  20  40  60 


Figure  3.  Arrival  Angles  vs  Depth  at  0.25  nmi  in  an  Isovelocity  Environment:  : 

(a)  Flat  and  (b)  Range-Dependent  Bathymetry 

Using  the  same  notation  as  before,  the  amplitude  of  the  Gaussian  ray  bundle  is  now  defined 
by 


'F  = 


Pv.<^\ 


■JlKi 


exp|-0.S[(z-z,)/c7,f  J, 


(22) 


where 


Pv.O  -''oPr,v.0^^v.O^O 


(23) 


depends  only  on  the  source,  and  where 
cr^  =^max(Az^,8A) 


(24) 
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is  an  effective  standard  deviation  or  half  beamwidth,  and  A  is  the  wavelength  at  the  field  point. 

The  factor  fiv,o  was  chosen  so  that  the  energy  within  a  geometric-acoustic  ray  tube  equals  the 
energy  within  a  Gaussian  ray  bundle.  Equation  (24)  was  found  by  accident,  as  will  be  described 
shortly. 

To  find  eigenrays,  classical  ray  models  typically  perform  an  iterative  search  on  test  rays  or 
interpolate  in  range  or  depth  between  test  rays.  The  solution  is  local.  If  equation  (17)  is  used,  a 
classical  eigenray  is  determined  from  a  test  ray  that  roughly  falls  within  Azv/2  of  the  field-point 
depth  at  the  field-point  range.  In  contrast,  ray  bundles  are  global,  affecting  all  depths  to  some 
degree.  It  is  assumed  that  the  Gaussian  eigenray  amplitude  is  formed  by  the  power  addition 

(25) 


of  bundles  of  the  same  type.  The  eigenray  source  angle  de^,  horizontal  slowness  pr,e,  vertical 
slowness  pz^e,  boundary  phase  shift  0e,  and  travel  time  Tg  are  obtained  from  the  weighted 


averages: 

V 

(26) 

Pr.e  =  X'Y."^^Pry  ^ 

V 

(27) 

p 

r  z,e  *■  e  /.  j  vr  z,v  •> 

V 

(28) 

V 

(29) 

'T  _  T’ 

(30) 

V 


To  obtain  a  more  accurate  expression  for  travel  time,  the  ray  bundle  travel  time  Tv  is 
extrapolated  to  the  field-point  depth  z  by 


Tz,v^T^  +Pz.v{^-^v) 


(31) 


in  equation  (30)  before  the  Gaussian  weighting  function  is  applied.  Justification  for  the 
extrapolation  in  equation  (31)  is  based  on  the  relationship 


(32) 


and  figure  4.  The  Gaussian  eigenray  amplitude 
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and  field-point  angle 


(33) 


(34) 


yield  the  complex  valued  eigenray  pressure 

Pe  =  A  exp(/®r^  +  zO  J .  -  (35) 

An  additional  complication  makes  it  necessary  to  compute  two  weighted  averages  for  each 
ray  type.  The  first  is  the  contribution  from  ray  bundles  with  Azy  >  0,  the  second  with  Azy  <  0. 

Ray  bundles  must  be  separated  in  order  to  simulate  the  reversals  in  test  ray  depth  versus  source 
angle.  These  reversals  are  caustics  that  may  be  formed  by  the  environment  or  by  numerical 
approximations.  In  the  isovelocity,  flat-bathymetry  case  shown  in  figure  3,  the  arrival  angle 
versus  depth  curve  is  monotonic  for  each  ray  type.  (It  was  convenient  to  display  arrival  angle 
versus  depth  instead  of  test  ray  depth  versus  source  angle.  The  two  are  equivalent  in  the 
isovelocity,  flat-bathymetry  case.).  In  contrast,  the  range-dependent  bathymetry  produces  many 
reversals. 

For  simplicity,  the  current  model  assumes  that  caustics  occur  at  turning  points.  Hence,  the 
ray-tracing  algorithm  adds  -7t/2  to  the  phase  shift  0y  at  each  upper  and  lower  vertex.  Asymptotic 
theory  dictates  that  the  eigenray  phase  shifts  0e  corresponding  to  the  Azy  >  0  and  Azy  <  0 
weighted  averages  be  separated  by  an  additional  zdl.  These  assumptions  are  not  valid  in  many 
environments,  and  the  coherent  predictions  are  subject  to  errors.  However,  since  the  primary 
purpose  of  the  model  is  to  simulate  high-frequency  reverberation  by  power  addition,  phase  errors 
are  a  secondary  concern. 


Figure  4.  CASS  Travel  Time  Interpolation 
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If  lSzvltS.dv,o  is  slowly  varying,  and  equation  (24)  is  simplified  to 
o,  =  Az,/2 


(36) 


then  classical  eigenrays  can  be  approximated  by  Gaussian  eigenrays  with  a  high  degree  of 
accuracy.  In  other  words,  Gaussian  eigenrays  can  be  made  as  accurate  as  classical  eigenrays. 
This  claim  is  suggested  by  figure  5,  which  shows  for  various  values  of  n. 

When  n  =  1,  the  ratio  reduces  to 


The  rectangular  curve  in  figure  5a, 

1  if  2',  -  Az,  /  2  <  z  <  Zj  +  Az,  /  2 

(38) 

0  otherwise 

provides  a  useful  comparison.  Equation  (23)  for  Pv,o  sets  the  area  under  the  Gaussian  curve  to  the 
area  under  the  rectangular  curve. 

According  to  equation  (37),  the  Gaussian  and  rectangular  curves  differ  by  a  factor  of 
-JiTk  «  0.80  or  -1  dB  at  z  =  zi.  The  factor  improves  to  1.01  or  0.06  dB  at  the  midpoint  depth  zj 
in  figure  5b.  Here,  the  contributions  fi'om  three  adjacent  bundles  are  added,  assuming  that  the 
increments  in  9v,o  and  Zyare  constant.  Figure  5c  shows  the  effects  of  doubling  the  number  of 
bundles  from  three  to  six  by  cutting  the  increment  in  Ovja  in  half  Although  the  overall  fit  is 
significantly  better,  fluctuations  in  the  vicinity  of  the  midpoint  depth  remain  at  0.06  dB.  Figure  5c 
suggests  that  the  Gaussian-induced  fluctuations  could  be  reduced  by  an  averaging  algorithm.  This 
is  borne  out  by  figure  6,  which  averages  two  results.  The  source  angles  of  the  second  are 
midpoints  of  the  first.  In  short,  energy  is  conserved  in  a  single  ray  bundle,  but  an  adequate 
description  of  the  acoustic  pressure  requires  the  superposition  of  several. 

The  standard  deviations  cjv  must  be  bounded  from  below  to  prevent  the  amplitudes  Py  from 
growing  too  large  in  the  vicinity  of  caustics.  This  bound  will  now  be  estimated  for  a  well-known 
test  case  developed  by  Pederson  and  Gordon  (reference  10)  to  investigate  extreme  downward 
refraction.  Figure  7  contains  their  sound  speed  profile  and  ray  diagram  for  a  0.9144-km 
(1000-yd)  source  depth.  The  inverse  sound  speed  squared  is  linear  in  depth,  so  one  can  express 
ray  trajectories  in  closed  form.  However,  the  results  below  do  not  use  this  feature  and  may  be 
contaminated  to  a  small  degree  by  numerical  artifacts.  It  should  also  be  noted  that  the  example 
does  not  represent  realistic  ocean  conditions  in  that  the  sound  speed  is  too  small  throughout  most 
of  the  water  column. 
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DEPTH  DEPTH 


RELATIVE  POWER 


RELATIVE  POWER 


Figure  5.  Power  Addition  of  Gaussian  Ray  Bundles  Normalized  by  the 
Geometric-Acoustic  Pressure  Squared 


SOUND  SPEED  ( km/s )  RANGE  (km) 


Figure  7.  Formation  of  a  Caustic  Under  Extreme  Downward  Refracting  Conditions: 

(a)  Sound  Speed  Profile  and  (b)  Ray  Diagram  for  a  0.9144-km  (1000-yd)  Source  Depth 

Figure  8  illustrates  the  effect  of  changing  the  lower  bound  of  a  on  propagation  loss  for  a 
0.73 15-km  (800-yd)  field-point  depth  in  the  vicinity  of  the  caustic  at  2.89  km.  Decibel  levels 
were  computed  by  coherent  eigenray  addition  for 

^max(Az^,A)  (a) 

^max(Az^.,8A)  (b) 

^  I  max(A2^,4;rA)  (c) .  (39) 

max(Az^  ,1 61)  (d) 

^max(A7^,32/l)  (e) 

In  figure  8a,  cTv  is  too  small.  Although  bounded,  the  level  is  too  high  in  the  illuminated  region  to 
the  left  of  the  caustic  and  falls  off  too  rapidly  in  the  shadow  zone  to  the  right  of  the  caustic.  The 
opposite  occurs  in  figure  8e,  where  <Jv  is  too  large;  now  the  decibel  level  is  too  small  at  the  caustic 
and  falls  off  too  slowly  in  the  shadow  zone.  Equation  (24)  is  currently  based  on  figure  8b  because 
it  currently  gives  the  best  overall  performance  near  the  caustic  for  both  coherent  and  power 
additions.  Evidence  to  be  provided  later  will  indicate  that  this  remarkably  simple  formula  is  far 
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2.83  2.84  2.85  2.86  2.87  2.88  2.89  2.90 
RRNGE  (kn) 


Figure  8.  Propagation  Loss  vs  Range  in  the  Vicinity  of  a  Caustic  for  a  0.9144-km  (1000-yd) 
Source  Depth,  0.7315-km  (800-yd)  Field-Point  Depth,  and  Various  Lower 
Bounds  of  a:  (a)  yj2,  (b)  dX,  (c)  2%X  (d)  8X  (e)  ItiX 


more  general  than  would  be  expected.  Note  that  figures  8c  and  8d  are  closer  to  the  generally 
accepted  answer  in  the  vicinity  of  2.86  km.  Since  this  range  is  dominated  by  the  ray  that  grazes 
the  ocean  surface,  it  is  likely  that  a  phase  error  was  produced  by  the  treatment  of  boundary 
reflections,  which  topic  is  treated  next. 

As  is  customary  in  typical  ray-tracing  models,  rays  are  attenuated  at  the  ocean  boundaries  by 
appropriate  reflection  coefficients.  The  cumulative  effect  is  contained  in  the  pressure  ratio 

Tv  =n^^.v  (40) 


and  phase  shift 

(41) 

where  /i,vand  denote  surface  and  bottom  reflection  coefficient  amplitudes,  respectively,  and 
and  ^6,  V  are  the  corresponding  phase  shifts.  Specific  examples  of  reflection  coefficients  appear 
in  later  sections. 

The  ocean  boundaries  impose  additional  requirements  on  the  numerical  computation  of  both 
classical  and  Gaussian  eigenrays.  Although  the  difficulties  discussed  below  stem  fi'om  tracing  rays 
to  a  constant  range,  it  is  felt  that  the  benefits  of  stepping  in  range  overshadow  any  disadvantages. 

First  consider  a  constant  sound  speed  environment.  The  range  r  and  field  point  to  source 
depth  z  -  zo  for  the  direct  path  are  related  by 

z-Zq  =  /'tan^o.  (42) 

Holding  r  constant,  the  differential 

dz  =  r  sec^  O^dO^  (43) 

provides  the  lower  bound 

Az,  >  0  (44) 

for  the  change  in  ray  depth  with  source  angle.  In  other  words,  if  the  angular  spacing  between 
tests  rays  is  too  large,  there  will  be  only  one  direct  path  at  the  range  of  interest.  It  wilt  not  be 
possible  to  compute  the  required  increments  Azv  and  v.o-  Realistic  ocean  conditions  create 
additional  complications. 

If  computer  run  time  is  not  a  concern,  and  the  sound  speed  and  bottom  are  well  behaved, 
one  can  make  ts.6  v.o-  sufficiently  small  that  several  rays  of  the  same  type  cross  the  range  of 
interest.  The  practical  approach  taken  here  uses 
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=  z  —  Az  /  2 


(45) 


=z^+^zjl 

instead  of  equations  (13)-(16)  to  set  the  bounding  depths  and  zf  where 

Az,  =zj-z,  (46) 

is  the  difference  between  the  bottom  and  surface  depths  at  range  r.  Since  adjacent  rays  of  the 
same  type  are  unavailable,  the  Gaussian  bundle  distributes  its  energy  normally  about  the  central 
ray  with  a  standard  deviation  of  half  the  water  column. 

It  is  often  possible  to  create  adjacent  rays  of  a  desired  type  by  unfolding  rays  of  an  adjacent 
type.  To  demonstrate  the  importance  of  these  virtual  rays,  consider  an  eigenray  that  is  incident  to 
a  field  point  slightly  above  the  ocean  bottom  as  shown  in  figure  9.  Circles  along  the  vertical 
dashed  line  represent  the  depths  at  which  a  fan  of  reflected  test  rays  cross  the  field  point  range  r, 
with  ray  v  being  the  closest  to  the  bottom.  Besides  denoting  range,  the  horizontal  axis  also  serves 
to  illustrate  the  relative  amplitude  of  a  Gaussian  ray  bundle  that  is  centered  at  z^  and  has  a 
standard  deviation  of  C7=  I  Azvl  /2.  The  unshaded  portion  represents  reflected  energy  that  does 
not  contribute  to  the  incident  eigenray.  Unless  unfolded,  power  due  to  the  shaded  portion  will  be 
lost  to  the  bottom,  and  the  eigenray's  amplitude  will  be  too  small.  By  projecting  the  bundle  along 
incident  rays,  it  becomes  clear  that  the  shaded  portion  corresponds  to  incident  energy. 


Figure  9.  Creation  of  Virtual  Rays  by  Unfolding  Classical  Rays  at  a  Boundary 
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Except  in  the  simplest  cases,  it  is  not  clear  how  the  virtual  ray  parameters  ffyja,  p'r,v,p'z,v, 
T'v,  'P'v  should  be  unfolded  from  their  real  counterparts.  To  avoid  the  cumbersome  task  of 
adding  or  removing  the  effects  of  boundary  reflection  coefficients  during  the  unfolding  process, 
the  current  procedure  sets 

^v,0  = 

P'r.vO^P*r,.> 


Pz^vO  Pz,v  y 

T'^T 

V  v> 

fi’ 


(47) 


where  the  asterisk  (*)  indicates  a  reflection  about  the  appropriate  boundary  and  p  is  the  ray  index 
nearest  to  v  such  that  ray  p  and  the  eigenray  are  of  the  same  type.  Errors  in  virtual  ray 
parameters  increase  with  distance  to  the  boundary.  However,  the  distance  from  the  virtual  ray  to 
the  field  point  also  increases,  so  that  the  effect  on  the  eigenray  decreases. 

For  example,  referring  once  again  to  figure  9,  one  sees  that  (p'r,v,p'z,^  are  the  direction 
numbers  of  the  dashed  line  for  ray  v.  The  virtual  ray  amplitude  and  phase  are  obtained  from  ray 
V-  1,  the  direct  path  ray  nearest  to  ray  v. 

A  few  other  special  cases  require  attention.  For  brevity,  the  discussion  is  limited  to  one 
such  case — ^nearfield  propagation.  For  the  amplitude  of  the  Gaussian  ray  bundle  Tv  to  approach 
spherical  spreading  as  the  field  point  (r,  z)  approaches  (0,  zo),  it  is  necessary  that  the  standard 
deviation 


CTy  oc 


(48) 


in  the  vicinity  of  the  source.  Unless  equation  (24)  is  adjusted,  equation  (44)  indicates  that  the  SA 
lower  bound  of  ctv  will  be  invoked  for  short  ranges,  such  that 


r< 


Az„ 


8;i 


A0 


v,0 


A0 


(49) 


t'.O 


The  current  adjustment  replaces  the  8A  in  equation  (24)  with 


mmJ 


nr 
kl80^ 


Ml. 


(50) 


One  degree  or  7t/180  radians  was  found  to  work  better  than  Ad^o  in  several  test  cases. 
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3.  MODEL  COMPARISONS 


This  section  contains  two  sets  of  propagation  loss  test  cases.  For  convenience,  the  Gaussian 
ray  bundle  based  computer  code  has  been  given  the  acronym  GRAB.  The  first  test  case  compares 
GRAB  results  with  those  of  the  Navy  standards  PE  V3.4  (reference  6)  and  ASTRAL  V4.2 
(reference  7)  at  25  Hz  and  10  kHz  for  a  classic  convergence  zone.  Then  GRAB  and  the  more 
academic  models  EFEPE  (reference  8)  and  OASES  (reference  9)  are  used  to  predict  propagation 
loss  at  1,000  Hz  for  a  shallow- water  environment.  If  GRAB  compares  well  at  these  lower 
frequencies,  it  should  also  perform  well  above  20  kHz,  the  high  fi-equencies  of  interest.  We  begin 
with  a  brief  overview  of  the  comparison  models. 

ASTRAL  V4.2,  a  current  Navy  standard  model  for  range-dependent  environments,  was 
designed  for  computational  speed,  frequencies  less  than  1  kHz  and  deep-water  scenarios. 
ASTRAL  bundles  energy  into  50-100  “shmodes"  and  provides  a  range-  and  depth-averaged 
intensity  solution.  (Shmodes  are  smooth  mode-like  functions  that  decay  exponentially  beneath 
their  turning  points.)  The  ASTRAL  treatment  of  bottom  interactions  uses  a  weighted  average  of 
a  reflected  and  a  refracted  path. 

The  parabolic  equation  (PE)  was  introduced  into  underwater  acoustics  in  the  early  1970s  by 
Hardin  and  Tapper!  (reference  11).  They  devised  an  efficient  numerical  scheme  based  on  Fourier 
transforms.  PE  techniques  have  grown  and  are  now  popular  for  solving  range-dependent 
propagation  problems.  The  Navy  standard  PE  V3.4  model,  a  direct  descendant  of  the  original 
Hardin  and  Tapper!  work,  uses  split-step  Fourier  transforms  in  depth  to  compute  the  coherent 
pressure  field.  Numerical  aliasing  and  phase  errors  limit  the  types  of  scenarios  that  can  be 
modeled.  In  particular,  shallow-water  environments  where  the  bottom  properties  represent  an 
abrupt  change  from  the  water  column  should  be  avoided.  As  is  typical  of  PE  models,  the 
computation  time  and  storage  requirements  increase  with  frequency. 

More  robust  PE  solutions,  such  as  Collins'  finite  element  EFEPE  model  (reference  12)  have 
been  designed  to  bypass  the  split-step  Fourier  depth  transform.  This  model  is  more  accurate  than 
the  Navy  standard  PE  V3.4  for  bottom  interacting  scenarios  and  can  treat  higher  angle  energy. 
However,  EFEPE  does  not  incorporate  volume  attenuation  into  the  water  column,  an  important 
consideration  at  high  frequencies.  The  stratified,  geoacoustic  bottom  description  must  include  a 
highly  absorbing  basement  layer  in  order  to  simulate  a  radiation  condition  without  introducing 
erroneous  reflections. 

OASES  (Ocean  Acoustic  and  Seismic  Exploration  Synthesis)  models  seismo-acoustic 
propagation  in  horizontally  stratified  waveguides  using  wave  number  integration  in  combination 
with  direct  global  matrix  techniques.  The  model  is  a  more  robust  version  of  SAFARI  (Seismo- 
Acoustic  Fast  field  Algorithm  for  Range  Independent  environments)  (reference  13).  OASES 
includes  shear  properties  in  a  layered,  geoacoustic  bottom  and  eliminates  aliasing  problems  by 
moving  the  integration  contour  onto  the  complex  plane,  i.e.,  introducing  artificial  attenuation. 

To  avoid  model/model  comparison  complications,  the  same  submodels  (surface  loss,  bottom 
loss,  etc.)  were  used  whenever  possible.  Since  the  bottom  descriptions  required  by  the  GRAB, 
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EFEPE  and  the  OASES  model  are  sufficiently  different,  it  was  necessary  to  invoke  an  OASES 
option  that  generates  a  plane  wave  coefficient  fi-om  bottom  layers  and  implement  these  results  into 
a  GRAB  bottom  loss  table.  Range-dependent  sound  speed  is  not  resolved  so  easily 
(reference  14).  OASES,  being  range  independent,  is  based  on  one  sound  speed  profile.  PE, 
ASTRAL,  and  EFEPE  transition  abruptly  from  one  profile  to  another  at  user-provided  ranges. 
Only  GRAB  models  the  sound  speed  continuously  with  range  as  described  in  section  2. 


3.1  CLASSIC  CONVERGENCE  ZONE 

GRAB  predictions  for  a  classic  convergence  zone  were  then  compared  with  those  of  the 
Navy  standard  models  ASTRAL  V4.2  and  PE  V3.4.  Both  standard  models  should  provide 
reasonable  results,  with  PE  being  the  more  accurate  of  the  two.  The  sound  speed  profile  and  ray 
trace  for  a  792.48-m  (2600-fl)  source  depth  and  13,716-m  (15,000-ft)  water  depth  appear  in 
figure  10.  GRAB,  PE,  and  ASTRAL  25-Hz  and  10-kHz  propagation  losses  are  contoured  in 
figures  11,  12,  and  13  respectively. 

To  highlight  the  formation  of  shadow  zones  and  caustics,  GRAB,  PE,  and  ASTRAL  were 
run  without  volume  attenuation  or  bottom  interacting  energy.  The  GRAB  prediction  used 
random  phase  eigenray  addition,  sampled  range  in  1.852-km  (1-nmi)  increments,  and  employed  a 
±20°  vertical  aperture  in  0. 1°  increments.  ASTRAL  also  used  a  1 .852-km  increment,  while  the 
full-field  PE  model  automatically  selected  sampling  parameters  to  avoid  aliasing. 
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Figure  10.  Classic  CZ  Environment:  (a)  Sound  Speed  Profile  and  (b)  Ray  Trace 
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Figure  11.  GRAB  Propagation  Loss  Contours  for  Classic  CZ  Environment: 

(a)  25  Hz  and  (b)  10  Hz 


IMll 
Reverse  Blank 


DEPTH  (KM) 


40 


160 


<75  78  81  84  87  90  93  96  99  102  105  108  >111  DB 

Figure  12.  PE  Propagation  Loss  Contours  for  Classic  CZ  Environment: 


(a)  25  Hz  and  (b)  10  kHz 
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Figure  13.  ASTRAL  Propagation  Loss  Contours  for  Classic  CZ  Environment: 

(a)  25  Hz  and  (b)  10  kHz 
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Figure  14  displays  GRAB,  PE,  and  ASTRAL  propagation  loss  versus  range  for  a  792.48-m 
(2600-ft)  receiver.  Although  GRAB  was  not  designed  for  frequencies  as  low  as  25  Hz,  the  model 
tracks  the  general  trend  very  well  and  maintains  the  energy  levels  in  the  tails  of  the  caustic.  Table 
1  compares  the  computation  times  of  each  model  using  a  50-MHz  386  IBM  personal  computer. 

A  large  percentage  of  the  GRAB  computation  time  is  a  result  of  input/output  operations. 

Table  1.  Comparison  of  Computation  Times 


Model 

Computation  Time 

25  Hz 

10  kHz 

GRAB 

4.5  minutes 

4.5  minutes 

PE 

6  seconds 

26  hours 

ASTRAL 

5  seconds 

5  seconds 

3.2  SHALLOW-WATER  ENVmONMENT  AT  1000  HZ 

The  next  set  of  test  cases  compares  GRAB,  EFEPE,  and  OASES  propagation-loss 
predictions  for  a  shallow- water  environment.  Due  to  the  high-frequency  computation  limitations 
of  EFEPE  and  OASES,  it  was  necessary  to  baseline  GRAB  at  1,000  Hz.  It  was  assumed  that  the 
essential  physics  at  1000  Hz  is  the  same  as  at  20  kHz,  and  that  differences  can  be  attributed  to 
surface,  volume,  and  boundary  losses  whose  accuracy  is  independent  of  the  propagation  model. 

For  the  model/model  comparisons  here,  the  environment  consists  of  the  sound  speed  profile 
in  figure  15,  a  perfectly  reflecting  sea  surface,  and  zero  volume  attenuation.  The  source  and 
receiver  depths  are  30.48  m  (100  ft).  Since  OASES  is  a  range-independent  model,  predictions 
involving  OASES  assume  a  flat  bottom  at  152.4  m  (500  ft)  of  all  sand  or  all  rock.  After  GRAB 
was  baselined  with  OASES  and  EFEPE  for  this  range-independent  environment,  GRAB  and 
EFEPE  predictions  were  compared  for  a  downslope  bathymetry.  Ray  diagrams  for  both  cases  are 
traced  in  figure  16. 


3.2.1  Sand  Bottom 

The  sand  bottom  parameters  include  a  1572-ni/sec  sediment  sound  speed,  a  1.268-gm/cm^ 
density,  and  an  attenuation  of  0.02  dB/A  or  0.03  dB  at  1,000  Hz.  These  yield  the  bottom  loss 
versus  grazing  angle  curves  in  the  upper  portion  of  figure  17.  GRAB  exercised  its  Rayleigh 
reflection  coefficient  option,  while  EFEPE  and  OASES  computed  equivalent  plane  wave 
reflection  coefficients  for  a  finite  and  halfspace  bottom,  respectively. 
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Figure  14.  GRAB,  PE,  and  ASTRAL  Propagation  Loss  Predictions  at  a  2600-ft  Receiver 
in  a  Classic  CZ  Environment:  (a)  25  Hz  and  (b)  10  kHz 
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Figure  1 5.  Shallow-  Water  Environment  Sound  Speed  Profile 

EFEPE  and  OASES  propagation  loss  predictions  are  compared  in  figure  18.  The  EFEPE 
prediction  used  a  0.1-m  depth  mesh  and  a  1.5-m  range  step.  Input  phase  velocities  for  the 
OASES  wave  number  integrand  spanned  1,000  to  8,000  m/sec,  with  the  latter  producing  a  1.1-m 
range  step.  Although  the  EFEPE  and  OASES  range  steps  are  different,  the  two  results  are  nearly 
identical. 

Figure  19  compares  the  coherent  GRAB  prediction  with  the  EFEPE  prediction  of  figure  18. 
The  GRAB  range  step  is  10~m  and  its  aperture  is  ±80°  with  test  rays  every  0.1°.  The  coherent 
GRAB  levels  differ  from  the  EFEPE  levels  (and  hence  the  OASES  levels)  by  less  than  1  dB  at 
most  ranges.  Note  the  exceptional  agreement  between  interference  patterns. 


3.2.2  Rock  Bottom 

The  rock  bottom  is  described  by  a  3750-m/sec  sound  speed,  a  0.02-dB//l  attenuation,  and 
2.5-gm/cm'’  sediment  density.  GRAB,  EFEPE,  and  OASES  bottom  loss  functions  for  a  rock 
bottom  are  shown  in  the  lower  portion  of  figure  17.  Note  that  the  sand  bottom  losses  for  grazing 
angles  beyond  20  are  highly  attenuated,  while  propagation  over  rock  is  strong  for  grazing  angles 
up  to  80°.  In  this  test  case,  an  80°  source  angle  interacted  with  the  bottom  50  times  in  2.6  km. 
Thus,  small  differences  in  the  bottom  description  become  significant. 
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Figure  16.  Shallow-Water  Environment  Ray  Trace:  (a)  Flat  152.4-m  (500 ft)  and 

(b)  Downslope  Bathymetry 
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Figure  1 7.  Bottom  Forward  Reflection  Loss  vs  Angle  at  1  kHz:  (a)  Sand  and  (b)  Rock 

Figure  20  illustrates  the  GRAB  random  and  coherent  propagation  loss  predictions  with  a 
1 . 1-m  range  step  and  a  ±80°  aperture  sampled  every  0. 1°.  Because  angles  up  to  80°  are  now 
important,  the  interference  pattern  is  more  complicated  than  it  was  for  the  sand  bottom.  To  make 
meaningful  comparisons,  the  intensity  is  averaged  over  a  0.1852-m  (0.1 -nmi)  sliding  range 
window.  Figure  21  displays  the  results  every  0.01852  m  (0.01  nmi).  As  in  the  sand  bottom  case, 
EFEPE  was  run  with  a  depth  mesh  of  0. 1  m  and  a  range  step  of  1 .5  m,  and  the  OASES  range  step 
for  the  phase  velocity  interval  1000  to  8000  m/sec  is  1 . 1  m.  Even  with  its  better  phase  modeling, 
EFEPE  results  have  phase  errors  due  to  high  propagation  angles.  Nevertheless,  the  differences 
among  the  range-averaged  GRAB,  EFEPE,  and  OASES  predictions  are  small. 
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3.2.3  Rock  to  Sand for  Variable  Bottom  Depth 

The  actual  scenario  that  prompted  this  series  of  tests  was  a  downslope  run  where  the 
bottom  properties  changed  from  rock  to  sand.  Figure  16  shows  that  the  slope  redirects  energy 
away  from  the  surface  duct.  OASES  is  no  longer  applicable  in  this  range-dependent  environment, 
and  GRAB  must  be  baselined  with  EFEPE  predictions. 

First  consider  the  effect  of  the  slope  with  an  all  sand  bottom.  Figure  22  shows  GRAB 
propagation  loss  contours  for  the  flat  152.4  m  and  downslope  bathymetries.  As  the  ray  trace 
figure  16  indicates,  the  overall  levels  at  the  30.48-m  receiver  are  lower  when  the  bottom  slopes 
downward. 

Note  that  the  coherent  propagation  losses  at  the  30.48-m  receiver  predicted  by  EFEPE  and 
GRAB  in  figure  23  transition  from  rapid  fluctuations  over  the  rock  bottom  to  slower  fluctuations 
over  sand.  Little  difference  is  seen  in  propagation  loss  levels  over  the  rock  (where  the  bottom  is 
still  relatively  flat),  but  significant  differences  are  seen  between  models  over  the  sloping  sand 
bottom  at  2.315  km.  Figure  24  compares  0.01852-m  range-averaged  EFEPE  and  GRAB 
predictions  with  the  random  phase  GRAB  prediction. 

Figure  25  shows  coherent  and  random  GRAB  predictions  at  21  kHz  without  volume 
attenuation.  The  1-kHz  and  21 -kHz  predictions  are  similar  because  the  bottom  loss  is  nearly  the 
same  for  the  two  frequencies,  and  neither  has  surface  or  volume  losses.  Thus,  the  comparison 
illustrates  how  energy  is  distributed  in  the  water  column.  The  higher  interference  rate  at  21  kHz 
gives  rise  to  rapid,  coherent  1 5-dB  fluctuations.  In  many  applications  the  random  addition 
provides  a  more  robust  measure  of  the  propagation  loss. 
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Figure  22.  GRAB  Propagation  Loss  Contours  at  1  kHz  for  a  Sandy  Bottom:  (a)  Flat 

and  (b)  Downslope  Bathymetry 
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Figure  23.  EFEPE  and  GRAB  Propagation  Loss  Predictions  at  1  kHz  for  a  Rock-to-Sand 

Downslope  Bathymetry 
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Figure  24.  EFEPE  and  GRAB  Range-Averaged  Propagation  Loss  Predictions  at  1  kHz 
for  a  Rock-to-Sand  Downslope  Bathymetry 
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Figure  25.  GRAB  Propagation  Loss  Predictions  at  21  kHz  for  a  Rock-to-Sand 

Downslope  Bathymetry 
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4.  SUMMARY  AND  CONCLUSIONS 


To  summarize,  a  high-frequency  propagation  loss  model  was  developed  to  investigate 
shallow  ocean  environments.  The  model  is  based  on  N  x  2D  ray  tracing,  Gaussian  ray  bundles, 
and  virtual  rays.  The  ray  tracing  determines  ray  trajectories,  inclination  angles  and  losses  due  to 
volume  attenuation,  and  the  ocean  boundaries  along  numerous  test  rays.  The  Gaussian  ray 
bundles  replace  the  classical  spreading  loss  of  geometrical  acoustics.  The  area  under  a  ray  bundle 
is  chosen  to  conserve  energy,  while  the  algorithm  for  its  standard  deviation  was  derived  by  fitting 
results  at  a  smooth  caustic.  Contributions  from  the  tails  of  bundles  that  extend  into  the  ocean 
boundaries  are  recovered  by  unfolding  test  rays  into  virtual  rays  at  the  appropriate  boundary. 

Ray  bundles  that  undergo  the  same  number  of  surface  and  bottom  reflections  are  combined 
into  acoustic  eigenrays.  An  eigenray  amplitude  is  simply  the  power  sum  of  its  ray  bundles. 
Weighted  sums  yield  the  remaining  parameters.  The  random  addition  propagation  loss  at  a  field 
point  is  the  power  addition  of  all  eigenrays  to  the  point  and  equals  the  power  addition  of  the 
contributing  ray  bundles.  However,  the  coherent  addition  propagation  loss  is  the  coherent 
addition  of  all  eigenrays  and  generally  does  not  equal  the  coherent  addition  of  the  ray  bundles. 

Although  rigorous  justification  is  not  currently  available,  the  Gaussian  ray  bundle  method 
appears  to  be  valid  over  a  larger  band  of  frequencies  than  originally  intended.  Its  applicability  is 
being  established  by  comparison  with  various  Navy  “standards"  and  at-sea  measurements.  Results 
for  the  classic  convergence  zone  and  shallow-water  test  cases  are  encouraging. 

Future  plans  include  resolving  some  of  the  more  questionable  aspects  of  the  model.  In 
addition  to  the  lack  of  mathematical  rigor,  phase  errors  occur  at  turning  points  and  caustics.  As 
with  many  ray-tracing  models,  the  selection  of  test  rays  is  left  to  the  user.  Too  few  may 
jeopardize  fidelity;  too  many  may  be  costly.  Current  experience  calls  for  vertical  angle  increments 
of  0.1  to  0.01°.  Away  from  caustics,  real  eigenrays  are  not  sensitive  to  the  lower  bound  of  the  ray 
bundle  standard  deviation.  Since  imaginary  eigenrays  are  formed  by  the  tails  of  Gaussian  bundles, 
and  these  tails  depend  on  the  standard  deviation,  one  should  expect  shadow  zone  propagation  to 
be  very  sensitive  to  the  lower. bound.  It  is  not  known  why  the  expression  for  the  lower  bound  is 
so  robust. 
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APPENDIX 

RAY-TRACING  EQUATIONS 


Ray-tracing  equations,  a  generalized  form  of  Snell's  Law,  may  be  obtained  from  the 
equation  of  characteristics  for  partial  differential  equations.  In  rectangular  (x,  y,  z)  coordinates. 
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where 


ds  =  -^[dxy  +{dyy  +{dzy 


(A-2) 


is  the  differential  arc  length  along  a  ray  path  and  the  sound  slowness  p  =  c'*  is  the  reciprocal  of 
sound  speed  c.  By  confining  the  rays  to  the  vertical  plane  y  =  0  and  setting  the  horizontal  range 
X  =  r,  equations  (A-1)  reduce  to 
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Multiplying  these  equations  by  p  and  applying  the  chain  rule  of  ordinary  differential  equations 
yields 
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or,  equivalently, 
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In  practice,  equations  (A-5)  require  numerical  integration.  Two  general  approaches  are  in 
use.  The  first  assumes  that  the  sound  speed  is  a  given  function  of  position  and  approximates  ray 
paths  by  evaluating  standard  integration  formulas.  The  second  approach  approximates  the  sound 
speed  in  sectors  so  that  equations  (A-5)  can  be  evaluated  in  closed  form.  The  former  approach 
tends  to  be  more  efficient  if  the  integration  step  size  is  smaller  than  the  distance  across  sectors. 
However,  because  it  is  difficult  to  estimate  a  suitable  step  size  for  the  realistic  conditions 
investigated  here,  the  latter  approach  was  selected  instead.  To  be  specific,  the  r  -z  plane  is 
divided  into  triangular  sectors,  in  each  of  which  the  gradient  of  slowness  squared  is 

=  constant.  (A-6) 

dz ) 


It  should  be  noted  that  the  ocean  sound  speed  depends  on  temperature,  salinity,  and 
pressure.  Since  the  variations  of  temperature  and  salinity  are  small  in  deep  isothermal  layers,  the 
dependence  of  sound  speed  will  be  nearly  linear  with  depth  because  of  the  increase  of  pressure 
with  depth.  In  this  case,  it  may  be  more  efficient  to  approximate  the  sound  speed  in  sectors  in 
which  the  sound  speed  gradient  rather  than  Vp^  is  constant.  The  approach  selected  here  offers 
certain  numerical  advantages  in  shallow  water,  such  as  working  with  parabolic  ray  segments 
instead  of  circular  arcs. 

Given  equation  (A-6),  one  can  integrate  equations  (A-5)  by  rotating  the  coordinate  system 
so  that  one  component  of  Vp^  vanishes.  The  implementation  used  here  integrates  equations 
(A-5)  directly  to 
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where  the  subscript  ( )c  refers  to  values  at  the  given  point  (r^,  Zc).  Upon  taking  square  roots,  one 
obtains 
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for  the  horizontal  and  vertical  components  of  slowness  along  the  ray  segment.  It  has  been 
assumed  that  the  range  is  strictly  increasing.  Also,  the  ±  sign  must  be  chosen  so  that  pdzids 
agrees  with  the  vertical  direction  of  the  ray  segment.  Reversals  in  direction  will  be  discussed 
shortly;  for  the  time  being,  assume  that  there  is  no  reversal  in  range  or  depth.  Then,  eliminating 
pids  from  equations  (A-8)  yields 
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Since  the  left-hand  side  of  equation  (A-9)  depends  explicitly  on  range,  while  the  right-hand  side 
depends  explicitly  on  depth  (range  and  depth  are  related  implicitly  by  the  ray  path),  equation 
(A-9)  may  be  integrated  from  (r^,  Zc)  to  (r^,  Zd)  to  yield 
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Rationalizing  the  numerators  of  equation  (A-10)  and  dividing  by  2  gives 
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where 


Ar  =  r^-r^, 


Az  =  z^-z^, 


(A-12) 


and  Ah  is  introduced  for  convenience.  Note  that  the  +  sign  has  been  removed  by  the 
rationalization  process. 

Two  final  auxiliary  formulas  are  required.  Setting  (r,  z)  to  (r^,  Zd)  in  equations  (A-7)  and 
using  the  differences  of  two  squares,  it  is  seen  that 
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Now,  suppose  that  rj  is  given  but  zj  is  currently  unknown.  Then,  (pdrlds)^  can  be  found 
from  the  first  of  equations  (A-8),  and  Ah  can  be  found  from  the  left-hand  side  of  equation  (A-1 1). 
Given  Ah,  (pdzlds)^  can  be  obtained  readily  from  the  second  of  equations  (A-13).  This  allows 
solution  of  equation  (A-1 1)  for 
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Similarly,  if  zj  is  given  but  id  is  unknown,  (pdzlds)d  can  be  found  from  the  second  of 
equations  (A-8),  and  Ah  can  be  found  from  the  middle  term  of  equation  (A-1 1).  Given 
Ah,{pdrlds)d  can  be  obtained  readily  from  the  first  of  equations  (A-13),  and  equation  (A-1 1) 
yields 


Ar  -  AH 


f 


d 


(A-15) 


This  completes  the  case  in  which  there  is  no  reversal  in  direction.  Now  suppose  that 
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To  prevent  the  physically  impossible 


the  ray  segment  is  terminated  at  the  vertex  depth  Zd  where 


(A-17) 
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That  is,  the -ray  is  horizontal  at  (r^,  z^).  In  the  next  case  to  be  considered,  the  ray  originates  at  a 
vertexing  depth.  Here,  the  vertical  direction  must  be  chosen  so  that 

(z-zj^>0.  (A-20) 

Finally,  if  both  (pdz/ds)c  and  vanish,  the  ray  is  a  horizontal  segment. 
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